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INTRODUCTION 


Tile  High-Frequency  Benchmark  is  a  computer  code  being  developed  to  employ  a 
sophisticated  ray  tracing  program  coupled  with  a  state-of-the-art  ionospheric  model  of 
the  user’s  choice.  The  merging  of  these  two  major  components  depends  critically  on 
the  development  of  a  spatial  smoothing  routine  that  accepts  arbitrary  ionospheric 
models  as  input.  Such  a  spatial  smoothing  routine  was  developed  by  Paul  (1991)  and 
implemented  in  the  Benchmark  (Ferguson  and  Shellman,  1991).  This  report  describes 
improvements  and  enhancements  to  that  computer  program. 

Most  models  of  high-frequency  (HF)  propagation  rely  on  synoptic  data  to  develop 
global  maps  of  key  propagation  parameters  such  as  the  maximum  useable  frequency 
(MUF).  One  failing  of  this  approach  is  that  large  errors  in  prediction  of  critical  propa¬ 
gation  parameters  occur  under  many  important  conditions.  This  is  especially  true  at 
high  latitudes.  In  part,  errors  in  the  models  may  be  attributed  to  the  fact  that  they  are 
based  almost  exclusively  on  hourly  averages  taken  at  mostly  middle  latitude  stations. 
Another  important  contribution  to  errors  is  the  lack  of  reliable  models  of  intermittent 
phenomena  such  as  sporadic-E.  The  most  widely  used  HF  prediction  codes  in  use  have 
no  provision  for  first  principles,  physics-based  calculations.  In  the  past,  first  principles 
codes  were  too  slow,  complex,  and  user-unfriendly  to  address  most  classes  of  problems 
of  interest  to  HF  communicators  and  model  developers.  Today,  modem  desktop  com¬ 
puters  rival  the  mainframes  of  the  past,  and  the  prospect  of  using  modem  first  princi¬ 
ples  codes  in  routine  propagation  problem  solving  is  a  reality.  Future  HF  prediction 
models  will  be  hybrids  of  rigorous  physics-based  propagation,  expert  systems  to  moni¬ 
tor  and  interpret  realtime  geophysical  and  ionospheric  sensors,  and  improved  man- 
machine  interfaces. 

The  goal  of  the  work  described  in  this  report  is  to  develop  techniques  that  permit 
timely  and  efficient  evaluation  of  different  propagation  models.  The  resulting  program 
is  to  be  as  sophisticated  as  possible  with  concerns  for  computer  run  time  to  be  sub¬ 
ordinate  to  the  use  of  realistic  models  of  die  environment.  Employed  with  new  data¬ 
bases  of  high-latitude  propagation  measurements,  this  model  will  be  used  to  evaluate 
and  improve  the  faster  running  models  and  point  to  deficiencies  in  the  ionospheric 
models.  In  addition,  the  development  of  new  graphical  displays  will  enhance  under¬ 
standing  of  the  propagation  environment. 

SPATIAL  SMOOTHING 

The  first  step  in  the  development  of  the  HF  Benchmark  was  to  develop  a 
generalized  spatial  smoothing  routine  that  would  accept  arbitrary  ionospheric  profiles. 

A  candidate  spatial  smoothing  algorithm  was  developed  by  Paul  (1991),  who  demon¬ 
strated  a  general  function  that  permitted  fitting  a  set  of  irregularly  spaced  discrete  data 
with  satisfying  results.  Ferguson  and  Shellman  (1991)  applied  this  algorithm  to  a  three- 
dimensional  array  of  ionospheric  profiles  of  electron  density  with  altitude,  distributed 
in  an  irregular  grid  of  latitude  and  longitude.  Modifications  have  been  made  to  that 
application.  For  completeness,  this  algorithm  will  be  reviewed. 
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The  basic  procedure  in  this  algorithm  is  to  define  an  ionospheric  parameter  at  a 
series  of  coordinates.  The  earlier  version  of  the  program  used  the  logarithm  of  the 
electron  density.  Currently,  the  program  is  using  the  logarithm  of  the  plasma 
frequency.  A  single  data  point  is  characterized  by  a  function,  / ,  where 


fi(x,y,2)  =  !+• 


a, 


1  + 


and  the  value  at  a  point  (x,y,z)  due  to  all  of  the  data  points  in  the  set  is  given  by 

p(x,y,2)  =  po+p,Y\fiix,y,2). 

The  quantities  r^i,  ryi,  and  are  called  the  control  distances  and  basically  determine 
the  range  of  influence  of  individual  data  points.  These  control  distances  should  be 
about  80%  of  the  separation  of  the  data  points  in  each  direction.  The  quantities  po  and 
Ps  are  scaling  parameters.  We  use  them  to  normalize  /  and  to  force  its  range  of 
values  to  be  between  0  and  1.  This  seems  to  improve  the  convergence  of  the  iterations 
used  to  find  the  values  of  the  unknown  coefficients,  Cj . 

Modifications  to  the  original  procedure  of  Paul  (1991)  were  required  to  consistently 
achieve  convergence  of  the  coefficients.  Given  that  there  is.  no  functional  form  to 
itnpose  on  the  iterations,  the  following  procedure  is  used.  The  minimum  value  of  the 
electron  density  data  is  determined  (po)  and  subtracted  from  all  of  the  values.  These 
data  are  then  normalized  using  the  peak  value  of  the  offset  data  (p,).  Equation  2  is 
used  to  compute  a  value  at  each  of  the  input  locations.  A  correction  term  is  defined 
equal  to  the  difference  between  the  computed  and  input  values  of  the  ionospheric 
parameter  at  the  input  locations.  If  the  largest  of  these  correction  terms  is  less  than  a 
specified  amount,  then  the  iterations  are  terminated.  Otherwise,  the  coefficients  are 
adjusted  by  some  fraction  and  the  calculations  repeated.  Currently,  we  are  using  a 
maximum  difference  between  these  calculations  and  the  input  data  of  0.05,  which 
amounts  to  5%  of  the  logarithm  of  the  peak  plasma  frequency  and  gives  reasonable 
results.  The  large  number  of  data  points  required  to  give  a  reasonable  representation 
of  the  ionospheric  variation  along  a  propagation  path  still  requires  a  large  number  of 
iterations.  Typically,  quite  a  few  iterations  are  required  to  resolve  individual  data 
points  and  that  the  convergence  is  slow  (many  iterations).  In  spite  of  this  slow  conver¬ 
gence,  using  a  larger  adjustment  factor  at  these  times  causes  the  iteration  process  to 
fail.  It  seems  that  we  must  always  use  a  small  fraction  to  adjust  the  coefficients 
(0.03125  in  the  current  configuration  of  the  program).  This  small  value  prevents  oscil¬ 
latory  behavior  of  the  iterative  corrections  but  requires  a  lot  of  iterations.  After  much 
experimentation,  we  found  that  the  initial  estimate  for  the  a/  should  be  some  fraction 
of  the  quantity  d/  -  1.  Currently  the  program  uses  0.6  for  this  fraction. 

Three  subroutines  are  used  to  implement  the  interpolation.  The  routine 
AKPJAODEL  takes  the  initial  path  parameters  and  sets  up  the  array  of  data  points  by 
calculating  the  locations  of  the  ionospheric  profiles.  The  routine  AKP_COEFF  finds  the 
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quantities  po  and  p,,  and  perfoims  the  iterations  to  determine  the  coefficients,  Ci.  The 
routine  AKP_XYTHPH  does  coordinate  conversion  between  the  magnetic  dipole  coordi¬ 
nate  system  used  in  the  ray  trace  routines  and  the  quasi-Cartesian  coordinate  system 
used  in  the  interpolation  routine.  This  routine  also  calculates  the  derivatives  of  the 
ionospheric  parameters  due  to  the  spatial  variation  of  the  ionosphere. 

IONOSPHERIC  MODEL 

The  ionospheric  model  currently  used  in  the  program  is  the  Utah  State  University 
(USU)  ionospheric  model,  a  parameterized  version  obtained  from  Anderson* .  This 
model  is  called  the  High  Latitude  Ionospheric  Specification  Model  (HLSM).  A  driver 
subroutine  named  USUJdODEL  has  been  developed  for  this  model.  This  routine  is 
called  by  AKP_MODEL  to  generate  ionospheric  profiles  at  points  500  km  apart  along 
the  propagation  path  (defined  by  the  transmitter  and  the  receiver)  and  at  points  500 
km  away  on  either  side  of  the  path  along  a  direction  perpendicular  to  the  propagation 
path.  The  previous  model  calculated  the  ionospheric  parameters  at  heights  determined 
by  the  HUSM  model  that  are  closer  together  at  the  bottom  of  the  ionosphere  than  at 
the  top.  The  heights  passed  to  the  spatial  smoothing  routine  were  evenly  spaced, 
having  been  interpolated  from  these  irregularly  spaced  values.  The  current  version  of 
the  program  now  uses  the  values  as  returned  from  the  ionospheric  model.  This  gives 
more  accurate  values  of  the  details  of  the  ionospheric  parameters  at  low  altitudes,  as 
originally  intended  for  the  HUSM. 

The  HUSM  is  used  because  it  is  considered  to  be  one  of  the  best  high-latitude 
models  available.  Unfortunately,  it  is  a  model  for  positive  ion  densities  in  the  altitude 
range  from  100  to  800  km.  To  get  electron  density,  charge  neutrality  is  invoked  and 
the  electron  density  is  set  equal  to  the  total  ion  density.  This  is  valid  only  where  the 
negative  ions  can  be  ignored.  The  model  tends  to  generate  (correctly)  large  ion 
densities  at  the  bottom  of  the  ionosphere  that  would  normally  be  neutralized  by  the 
presence  of  negative  ions  as  well  as  electrons  at  those  altitudes.  For  the  moment,  an 
adjustment  to  the  model  has  been  made  to  force  the  electron  density  to  decrease  with 
altitude  below  the  ionosphere.  The  user  selects  a  minimum  height  at  which  the 
ionospheric  density  is  computed.  Two  additional  heights,  separated  by  20  km,  are 
added  below  this  height  that  force  a  decade  decrease  in  the  plasma  frequency  at  each 
lower  (additional)  altitude.  Under  nighttime  conditions,  these  values  may  cause  a  slight 
kink  in  the  electron  density  profile.  The  adjustment  is  required  to  provide  a  gradual 
entry  into  the  ionosphere  for  the  ray  tracing  routine.  The  unrealistically  high  values  of 
electron  densities  at  the  bottom  of  the  ionosphere  were  causing  numerical  problems  in 
the  ray  tracing.  We  plan  to  obtain  and  implement  an  electron  density  model  either  as  a 
modification  of  the  USU  model  or  as  a  replacement.  A  new  ionospheric  model  based 
on  the  Parameterized  Realtime  Ionospheric  Specification  Model  (PRISM)**  will  be  im¬ 
plemented.  The  new  model,  the  Parameterized  Ionospheric  Model  (PIM),  differs  from 
PRISM  in  not  being  realtime.  However,  it  is  a  global  model  and  will  give  the  Bench¬ 
mark  a  global  ionospheric  model  suitable  for  a  variety  of  propagation  studies. 

*  Private  commuoication  with  D.  Aadenon,  1991. 

**  Private  communicatJoii  with  Daniels,  1992. 
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RAY  TRACING 


As  in  Ferguson  and  Shellman  (1991),  the  HF  Benchmark  continues  to  use  the  Jones 
and  Stephenson  (1975)  ray  tracing  model.  This  is  a  versatile  program  with  full  allow¬ 
ance  for  externally  specified  models  of  the  electron  density,  collision  frequency,  and 
geomagnetic  field.  Every  effort  has  been  made  to  retain  this  flexibility  in  its  implemen¬ 
tation  in  the  Benchmark.  In  particular,  the  input  requirements  are  nearly  identical  to 
the  original  Jones  and  Stephenson  version,  with  the  spatial  smoothing  algorithm  being 
specified  as  an  ionospheric  model  named  “AKPMODEL”  When  the  program  sees  • 

'*AKPMODEL”  as  the  ionospheric  model,  the  next  input  record  is  read  as  a  case  identi¬ 
fication,  for  example  CASE-ID,  and  forms  the  root  of  the  file  name  for  all  subsequent  • 

outputs.  For  example,  the  parameters  of  the  ray  paths  are  stored  in  a  file  of  the  name 
CASE-ID. TMP,”  and  the  file  named  “CASE-ID.002”  contains  an  information  log.  In 
addition  to  the  original  program  inputs,  additional  parameters  have  been  introduced. 

These  parameters  are  summarized  in  table  1.  If  the  value  of  parameter  100  is  9,  the 
program  looks  for  a  file  with  the  name  of  the  case  and  the  extension  AKP.  The  coeffi¬ 
cients  and  other  parameters  are  read  from  this  file.  This  is  illustrated  in  the  sample 
input  file  shown  below. 


Table  1.  Ray  trace  inputs  for  AKP  model. 


Index 

Value 

Usage 

100 

1 

This  value  used  to  signal  that  the  at  need  to  be  computed 
(9  indicates  that  the  o;  have  been  computed  previously). 

111 

1 

Number  of  the  month  (Jan  is  1). 

112 

31 

Day  of  the  month. 

113 

1991 

Year. 

114 

1200 

UT  in  HHMM  format. 

115 

100 

10-cm  flux. 

116 

1 

k. 

117 

1 

Sign  of  By  (1  for  -1  for  “-”), 

118 

120 

Minimum  USU_MODEL  height  to  use  (km). 

119 

400 

Maximum  VSUJAODEL  height  to  use  (km). 

120 

0.8 

Horizontal  control  distance  factor. 

121 

0.8 

Vertical  control  distance  factor. 

122 

0.6 

Initial  scale  factor  for  a^. 

123 

6.0 

Range  of  influence  factor. 

124 

0.03125 

Iteration  adjustment  factor. 

125 

0.05 

Maximum  error  to  terminate  iterations. 

126 

501 

Maximum  number  of  iterations. 

130 

Plot  of  ionospheric  parameters;  =0  no  output;  =1  logarithm  of 
electron  density;  =2  plasma  frequency;  =3  X=l,  X+Y=l,  X-Y=l. 

131 

2 

Plasma  frequency. 

132 

200 

Height  at  which  i.  do  the  horizontal  contours  of  the  ionospheric 
parameter. 
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Minor  procedural  modifications  have  been  made  to  facilitate  running  this  program 
on  a  PC  but  the  FORTRAN  used  in  the  program  is  not  specific  to  that  platform.  In 
fact,  an  identical  capability  exists  for  a  VAX. 

SAMPLE  PROBLEM 

A  sample  problem  is  presented  here  to  illustrate  the  unification  of  the  individual 
models.  The  path  is  defined  in  geographic  coordinates  from  68.5®N,  32.5“E  to  SS.S^N, 
ilT.l^W.  The  geographic  bearing  angle  is  339.3“,  and  the  path  length  is  6000  km.  The 
ionospheric  profiles  are  generated  for  23  December  at  0215  UT.  The  10-cm  flux  is  70, 
the  kp  is  2,  and  the  direction  of  the  sun’s  magnetic  field,  By,  is  positive.  The  input 
parameters  for  the  program  are  shown  in  figure  1.  Contour  plots  of  the  plasma  fre¬ 
quency  in  the  vertical  plane  along  the  path  are  shown  in  the  top  panel  of  figure  2,  and 
contours  of  the  plasma  frequency  in  the  horizontal  plane  at  an  altitude  of  200  km  are 
shown  in  the  bottom  panel.  This  display  is  generated  by  an  auxiliary  program  named 
PLOTPATH.  This  program  uses  the  same  input  as  the  Benchmark  and  calculates  the 
selected  ionospheric  parameter  using  the  file  containing  the  calculated  coefficients.  It  is 
clear  that  the  ionosphere  is  changing  in  all  three  dimensions.  In  particular,  we  note  a 
broad  maximum  near  2000  km  followed  by  a  constriction  of  the  contour  lines  near 
3500  km. 

Ray  paths  for  the  ordinary  ray  at  3  MHz  are  shown  in  figure  3.  The  initial  eleva¬ 
tion  angles  are  20“  and  50“.  The  top  panel  shows  the  projection  of  the  ray  paths  onto 
the  vertical  plane  passing  through  the  transmitter  and  the  receiver.  The  ray  with  a 
small  circle  on  its  end  (elevation  50“)  is  a  ray  that  penetrates  the  ionosphere.  The 
bottom  panel  of  figure  2  shows  the  projection  of  the  ray  paths  onto  the  ground  with 
locations  above  the  horizontal  axis  being  to  the  west  of  the  propagation  path.  We  see 
that  the  off-path  variation  of  the  ionosphere  causes  some  deflection  of  the  ray  paths. 
Ray  paths  for  the  extraordinary  ray  at  3  MHz  are  shown  in  figure  4.  As  in  figure  3, 
the  elevation  angles  are  20“  and  50“.  However,  the  ray  at  50“  does  not  penetrate  but 
does  show  a  strong  deflection  from  the  direct  path.  Figures  5  and  6  show  ray  traces 
for  this  path  at  4  MHz  at  elevation  angles  of  0*  and  30“.  We  see  that  the  ray  at  0“ 
becomes  a  grazing  ray  and  never  reaches  the  earth’s  surface.  The  ray  at  30“  shows  a 
very  nonuniform  pattern  of  earth  and  ionospheric  reflections,  confirming  the  non- 
uniform  ionosphere  indicated  in  figure  2. 

IONOSPHERIC  DISPLAYS 

Two  ionospheric  displays  have  been  developed.  One  of  the  displays  is  a  global  view 
showing  contours  of  the  ionospheric  profile  in  a  fixed  meridian,  and  the  other  is  a 
view  along  short  paths  (a  subset  of  the  global  view).  In  each  display,  contours  of  the 
electron  density  are  shown,  taken  directly  from  the  ionospheric  model.  Figure  7  shows 
an  example  of  the  global  display.  The  map  projection  is  orthographic,  with  the  obser¬ 
vation  point  in  the  center  of  the  display.  The  location  of  this  point  is  specified  by  the 
user.  Superimposed  on  the  circle  representing  the  earth  is  a  series  of  lines  indicating 
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land  masses.  There  is  also  a  series  of  lines  ranging  from  solid,  through  dashes  of 
various  lengths,  to  dotted.  These  lines  indicate  the  position  of  solar  zenith  angles 
ranging  from  90"  (solid)  to  99"  (dotted)  to  show  where  the  day-night  transition  is 
located.  If  it  is  visible  from  the  observation  side  of  the  globe,  the  sub-solar  point  is 
also  shown.  The  vertical  scale  for  points  above  the  earth’s  surface  is  specified  by  the 
user  by  setting  the  number  of  inches  between  the  earth’s  surface  and  the  top  of  the 
model  ionosphere.  Because  the  model  ionosphere  is  strictly  valid  at  high  latitudes,  the 
strange  features  shown  at  middle  and  equatorial  latitudes  must  be  ignored.  Neverthe¬ 
less,  the  auroral  bulges  are  visible,  and  the  differences  between  the  daytime  and  night¬ 
time  ionospheres  are  evident.  The  second  display,  shown  in  figure  8,  is  a  modification 
of  the  first.  The  primary  difference  from  the  global  display  is  the  shorter  path  length, 
which  allows  scaling  of  the  vertical  dimension  closer  to  the  horizontal  dimension.  The 
user  specifies  all  other  parameters. 


PROBLEMS  AND  FUTURE  PLANS 


Unfortunately,  the  main  problem  noted  for  the  first  program  remains;  namely,  the 
iteration  process  required  to  obtain  the  coefficients  for  the  interpolation  model  is 
excessive  in  computer  time,  and  the  interpolation  is  quite  slow.  A  faster  method  is 
highly  desirable.  ’The  current  interpolation  model  still  requires  too  much  manual  adjust¬ 
ment  to  be  routinely  used.  The  ionospheric  model  needs  to  be  changed  to  give  more 
realistic  electron  densities  in  the  lower  altitudes.  ’This  problem  will  be  addressed  by  the 
acquisition  of  the  newer  ionospheric  model  PIM.  This  latter  effort  is  very  important  if 
the  Benchmark  program  is  to  properly  account  for  absorption.  The  ray  trace  model  has 
some  difficulty  with  the  radical  changes  introduced  by  the  ionospheric  model.  Modifi¬ 
cations  have  been  made  to  the  ray  tracing,  but  additional  refinement  of  its  procedures 
are  probably  necessary. 
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Figure  1.  Input  parameters  for  the  sample  (transpolar)  path 
(continued). 
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Figure  3.  Plots  of  Uie  ray  paths  for  the  ordinary  ray  at  3  MHz  in  the  vertical  plane  passing  through  the 
transmitter  and  the  receiver  (top  panel)  and  in  the  horizontal  plane  (bottom  panel).  The  elevation  angles 
are  20*  and  50". 


Figure  4.  Plots  of  the  ray  paths  for  the  extraordinary  ray  at  3  MHz  in  the  vertical  plane  passing  through  the 
transmitter  and  the  receiver  (top  panel)  and  in  the  horizontal  plane  (bottom  panel).  The  elevation  angles  are 


Figure  5.  Plots  of  the  ray  paths  for  the  ordinary  ray  at  4  MHz  in  the  vertical  plane  passing  through  the  transmitter 
and  the  receiver  (top  panel)  and  in  the  horizontal  plane  (bottom  panel).  The  elevation  angles  are  0®  and  30®. 


Figure  6.  Plots  of  the  ray  paths  for  the  extraordinary  ray  at  4  MHz  in  the  vertical  plane  passing  through  the  trans 
mitter  and  the  receiver  (top  panel)  and  in  the  horizontal  plane  (bottom  panel).  The  elevation  angles  are  0®  and  30' 
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a  transmitter  and  a  receiver. 
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